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Abstract 



The reconstruction of the exchange-correlation potential from accurate ab initio electron densi- 
ties can provide insights into the limitations of the currently available approximate functionals 
and provide guidance for devising improved approximations for density- functional theory (DFT). 
For open-shell systems, the spin density is introduced as an additional fundamental variable in 
Spin-DFT. Here, we consider the reconstruction of the corresponding unrestricted Kohn-Sham 
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potentials from accurate ab initio spin densities. In particular, we investigate whether it is pos- 
sible to reconstruct the spin exchange-correlation potential, which determines the spin density in 

>: 

spin- unrestricted Kohn-Sham-DFT, despite the numerical difficulties inherent to the optimization 

O 

of potentials with finite orbital basis sets. We find that the recently developed scheme for unam- 
biguously singling out an optimal optimized potential [J. Chem. Phys., 135, 244102 (2011)] can 
provide such spin potentials accurately. This is demonstrated for two test cases, the lithium atom 
and the dioxygen molecule, and target (spin) densities from Full-CI and CASSCF calculations, 
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respectively. 
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1 Introduction 



Density-functional theory (DFT) within the Kohn-Sham (KS) framework [UI2] represents 
one of the most frequently applied quantum-chemical methods for electronic structure calcu- 
lations and for the determination of molecular properties. Its success relies on the accuracy 
of existing approximations to the exchange-correlation energy functional E xc [p] and to the 
exchange-correlation potential v xc [p] = 5E xc [p]/5p, i.e., the functional derivative of E xc [p] 
with respect to the electron density p(r) [3]. However, for open-shell systems, in particular 
for transition metal complexes, the existing approximations have a number of severe short- 
comings [1H6] , for instance for the prediction of the energy differences between different spin 
states [7HTT] and of spin-density distributions [T2HT5] . 

While the universal functionals E xc [p] and v xc [p] are unknown, there exists a numerical 
recipe for obtaining the exact exchange-correlation potential f xc [Po] corresponding to the 
ground-state electron density po of arbitrary atomic and molecular systems. First, this 
ground-state electron density po can be calculated accurately — and in principle exactly - 
using wave-function based ab initio calculations. Second, the Kohn-Sham potential v s [po\ 
that yields the density po i n a noninteracting system can be reconstructed. Finally, by 
subtracting the known nuclear and Coulomb potentials from this reconstructed potential, 
the exchange-correlation potential f xc [po] can be obtained. Such reconstructed ground-state 
exchange-correlation potentials can provide guidance for the construction of approximate 
exchange-correlation potentials [16-19J and energy functionals [20H22]. 

The key step in the above recipe is the reconstruction of the potential v s [po] from the target 
density po- This step corresponds to an inverse problem in quantum chemistry [23], i.e., 
the potential is sought which generates a given target density in a noninteracting reference 
system. This potential reconstruction is also essential for quantum-chemical subsystem and 
embedding methods (for a review, see Ref. [21]), in which it can be used to avoid the need 
for approximating the nonadditive kinetic-energy [25H2H], or for developing better approxi- 
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mations for this part of the embedding potential [3"0H3"3"] . 

The inverse problem of reconstructing the noninteracting local potential yielding a given 
target density is equivalent to evaluating the functional derivative of the noninteracting 
kinetic-energy functional T s [p] [3^135]. which is an implicit functional of the electron density. 
For evaluating such functional derivatives of implicit functionals, the optimized effective po- 
tential (OEP) method can be employed, which tackles the inverse problem by minimizing the 
implicit functional with respect to the local potential [36H38] . possibly subject to additional 
constraints [39J. Thus, the reconstruction of local potentials is a special case of the more 
general problem of evaluating the functional derivative of implicit density functionals. 

For applying OEP methods to many-electron systems and large molecules, both the orbitals 
and the local potential have to be expanded in finite basis sets. The introduction of a finite 
orbital basis set, however, turns the OEP method into an ill-posed problem and the solution 
becomes non-unique [30] • The ill-posed nature is common to many inverse problems and 
makes the inverse mapping from an electron density to a local potential unstable and sensitive 
to optimization parameters. Furthermore, these drawbacks result in unphysical potentials, 
which can contain large oscillations affecting orbital energies and derived properties [lOjUT]. 

To allow for a routine application of the OEP method in quantum chemistry, approaches to 
regularize the OEP solutions have been proposed. One approach developed by HeBelmann 
etal. is based on explicitly constructing an orbital basis set that is balanced with respect 
to the basis set employed for expanding the potential [32], whereas the orbital basis set is 
balanced implicitly in the approach of Kollmar and Filatov [13"|H3] . However, these methods 
require very large orbital basis sets, which hampers their application to larger molecular sys- 
tems. A different approach was developed by Yang and co-workers [3UI35], who introduced a 
regularization parameter in the energy functional to make the resulting optimized potentials 
are as smooth as possible. Recently, an approach which yields unambiguous potentials for 
any combination of orbital and potential basis sets and that provides high-quality poten- 



3 



tials already with small orbital basis sets was suggested by one of us [46J. It is based on 
the condition that the optimal reconstructed potential should yield the target density when 
extending the orbital basis set. 

So far, OEP methods for reconstruction the exchange-correlation potential from accurate 
ab initio densities have mainly been applied to closed-shell systems, i.e., to the total elec- 
tron density as target only (for exceptions, see, Refs. [311117]). For open-shell systems, one 
commonly employs an unrestricted KS-DFT formalism [TJ|l8j|l9], in which the spin den- 
sity Q(r) = p a (r) — p /3 (r) is used as an additional fundamental variable. This leads to 
separate KS equations for a- and /3-electrons containing different exchange-correlation po- 
tentials v° c [p,Q] = 5E xc [p,Q]/5p a and v^ c [p,Q] = 6E xc [p,Q]/5pp, i.e., the a- and /3-electron 
densities p a and // , respectively, are determined separately. In such a formalism, the exact 
spin-resolved exchange-correlation functional would yield — in addition to the exact total 
electron density — also the exact spin density [49H5Tj . While the total electron density is 
determined by the total exchange-correlation potential v^[p,Q) = \{v xc [p,Q] + v^ c [p, Q)), 
the spin density in unrestricted KS-DFT is determined by the spin exchange-correlation 
potential v x ^. m [p, Q] = ^(v xc [p,Q] — v^ c [p,Q]) [49J. Thus, for improving the spin-density de- 
pendance of approximate exchange-correlation potentials, it would be desirable to be able 
to reconstruct this spin exchange-correlation potential ^"[poj Qo] from accurate ab initio 
total and spin densities. 

When using finite orbital basis sets, such a reconstruction of v^ in [p , Q ] is particularly 
challenging, as it will be more sensitive to numerical errors than for the total or the individual 
a- and /3-electron potentials. Here, we extent the unambiguous potential reconstruction 
developed in Ref. (IB] to the spin-unrestricted cases and apply it to the reconstruction of the 
spin exchange-correlation potential from accurate ab initio (spin) densities. This requires 
some numerical enhancements of our implementation prompted by the use of Gaussian- 
type orbitals (GTOs) for expanding the orbitals in such wave-function based calculations. 
Moreover, the quality of the reconstructed spin potentials needs to be assessed carefully, if 
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possible by comparison to reference potentials reconstructed in fully numerical calculations. 

This work is organized as follows. Section |2] briefly reviews the potential reconstruction 
algorithm and outlines its extension to the spin-unrestricted case. In Section [3j the compu- 
tational methodology and extensions of our implementation are described. Subsequently, we 
study the reconstructed exchange-correlation potentials for two test cases, the lithium atom 
and the O2 molecule, in Section HI Here, target densities from both KS-DFT calculations 
and from accurate wave-function based ab initio calculations (Full-CI and CASSCF) are 
employed. Finally, our conclusions are summarized in Section 



2 Theoretical Background 



2.1 Determining Optimized Unrestricted Kohn— Sham Potentials 



Within spin-unrestricted KS-DFT, the wavefunction of the KS reference system is given by a 
(spin-unrestricted) iV-electron Slater determinant, which is constructed from N = N a + N 13 
orthonormal one-particle functions {<f)f(r)a(s), a = a, /?}. The corresponding spatial orbitals 
</>f (r) can then by determined by solving two separate sets of one-electron equations [19], 



(r) = eJW(r) 



and 



■-A + «f(r) 



and the a- and /3-electron densities are given by 

N a TV 3 



E 



and 



/(r) 



E 



(2) 



Here, we consider the inverse problem of determining the spin- resolved KS potential [i.e., 
the local potentials vf(r) and i>j (r)] from given a- and /3-electron target densities, Pq (r) and 
Pq(t), that is, we require 



p a (v) = pUr) and /(r) = pg(r) 



(3) 
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An alternative way of expressing this problem is to consider the total and spin densities [49J , 
p (r) = p a (r)+p /3 (r) and Q(r) = p Q (r) - /(r) (4) 

as target, and to regard the total and spin potentials [49J, 

*T M = ~ «(r) + «f (r)) and ^ in (r) = \ «(r) - «f(r)) (5) 

as the quantities that are sought. Here, f* ot (r) is the potential determining the total electron 
density, whereas w^ pm (r) determines the spin density. Therefore, this representation will be 
particularly useful for understanding the dependence the exchange-correlation potential on 
the spin density Q(r) and to identify the reason for the failure of approximate exchange- 
correlation functionals to describe the spin density correctly in some cases [15] . 

In principle, any method applicable for reconstructing the KS-potential in closed-shell sys- 
tems could be adapted to the spin-unrestricted case by applying it separately to the a- and 
/3-electron densities. However, already the closed-shell cases poses many numerical difficul- 
ties, and achieving uniform accuracy for the a- and /3-spin potentials, as it is required for 
obtaining t^ pin (r) accurately, turns out to be a challenging task. 

The conceptually simplest approach for determining the local potential yielding a given target 
density is to represent the potential numerically on a grid and to determine it iteratively. 
Several methods working along these lines have been developed over the past decades [HH 
[52H55] . Generally, these methods calculate the density from some trial potential and then 
update the potential by comparing the density to the target density. If the density is too 
large at a grid point, the potential is made more repulsive at this point. Conversely, if 
the density is too small, the potential is made more attractive. This process is repeated 
iteratively until the target density is obtained. Different numerical potential reconstruction 
methods differ in the way in which the potential is updated in each iteration [TB IIo'4'llo'o*] . The 
only exception is the method of Zhao-Morrison-Parr (ZMP) [53J, which uses a conceptually 
different approach. 
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However, such numerical methods also require that the KS equations [Eq. ([I])] are solved 
numerically on the same grid. Therefore, their application has mainly been limited to (closed- 
shell) atoms and, in some cases, (closed-shell) diatomic molecules [561157]. Here, we will 
employ such a fully numerical scheme to obtain accurate reference potentials for atoms. For 
determining optimized KS potentials in general molecular systems, both the orbitals and the 
potential are usually expanded in a basis set |371|39]. However, as will be discussed below, 
with finite orbital basis sets the potential reconstruction turns into an ill-posed problem, in 
which the resulting potential is not unique [HJJUTIEB] • 

To overcome the resulting numerical difficulties, we will apply the recently developed un- 
ambiguous optimization method and generalize it to unrestricted KS potentials. This 
scheme is based on a two-step procedure, in which one first determines a non-unique poten- 
tial using a direct optimization in a finite basis set. Subsequently, an unambiguous optimized 
potential is singled out by means of a suitable criterion. 

2.2 Direct Optimization of Unrestricted Kohn— Sham Potentials 

In the first step, two non-unique local potentials vf(r) and t>f (r) yielding the target a- 
and /3-electron densities pg (r) and Po( r )' respectively, in a given finite basis set have to be 
determined. To this end, we apply the direct optimization method by Wu and Yang [39J and 
extend it to the spin- unrestricted case. 

The KS kinetic energy for any pair of a- and /3-electron densities Pq,Pq is defined as [UG3] 
T s [plp%] = min mf |*.) = \t s [2 P «] + \t s [2 P 1] (6) 

where T = —A/2 is the kinetic-energy operator and where the minimization includes all 
wavefunctions ^ s corresponding to a spin-unrestricted iV-electron Slater determinant with 
a- and /3-electron densities p a (r) and p /3 (r). Hence, T s [pq,Pq] corresponds to the minimum 
kinetic energy of an unrestricted KS wave function \l/ s under the constraint that its a- and 
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/3-electron densities equal the target densities 

The constraint minimization problem of Eq. fl6]) can be reduced to two separate problems 
for 2p a and 2p@, respectively, which results in two Lagrangian functionals, jy Q [p a (r)] and 
K^[/r(r)], subject to the constraints of Eq. ([3]) with two corresponding Lagrangian multiplier 



functions, vf(r) and t>f (r), 

JV CT 

W a K\ = j>?|T|#> + I <(r)(p ff (r) - pt(r)) d 3 r for a = a,0. (7) 
i J 

Following Ref. [391 the local potentials which yield the target a- and /3-electron densities 
can now be determined by the unconstrained maximization of W /<J [p°"(r)] with respect to the 
local potential uj(r) for each electron spin a. 

To perform this maximization, the local potential is expanded in a finite basis set as [38|l39] , 

<(r) = i; cxt (r) + ^couiM(r) + v (r) +J2 b t9t(r), (8) 

t 

where t> nU c( r ) is the nuclear potential, vcoui( r ) is the Coulomb potential of the target density 
Po — Po + Po> an d v o( T ) represents an initial guess for the exchange-correlation potential, 
while the remainder is expressed as a linear combination of a finite set of basis functions 
{<7t(r)} with coefficients {b^}. For fixed t> ex t(r) arid vq(t), the unconstrained maximization 
of turns into an extremum problem with respect to the expansion coefficients 

for each electron spin. The first and second derivatives of ^""[/^(r)] with respect to 
can be calculated analytically and one obtains [39] the following expression for the gradient, 
BW a f 

-gyT = J 9t (r)(p°(r)-pUr))d 3 r (9) 
and the Hessian, 

st i a % a 

for each electron spin a. Note that Eqs. ([9]) and ffTOl are simplified for the case of real- valued 
orbitals here. With the gradient and Hessian available, the maximization can be performed 
using a standard Newton-Raphson optimization. 
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If a finite basis sets is employed for representing the KS orbitals, the potential reconstruction 
turns into an ill-posed problem and the optimized potentials resulting from the Wu-Yang 
direct optimization as described here are not unique [HUSHES]. This can be seen [16] by 
considering a change in the local a- or /3-electron potential Av J(r) = v°(r) — t>J (r), where 
Vg (r) is the potential obtained from the direct optimization, generating the orbitals {<P1 } 
and {4>a}- To first order, this change Av J (r) induces a response in the density, 



occ unocc 



= 2 E E ^"' A 5'5 ) C(r)C(r), (11) 

i a e i e a 

with Af CT (r) = ~Y^ t Ab%gt(r) and A6f = b a t — b1 . Here, one notices that any change in the 
potential AfJ(r) will leave the electron density unchanged if (0£|Ai£(r)|0f) = 0. Hence, if 
the orbital basis is not flexible enough, the a- and /9-electron densities are not affected by 
certain changes, e.g., oscillations, in the respective potentials. Linear combinations of basis 
functions <7t(r) for which the condition (0^|A-uJ (r)|0°") = holds are obtained by inserting 
the basis set expansion for the potential and performing a singular value decomposition 
(SVD) of the matrix B^ t = (4> a a (pi\gt) / - <0> which leads to 

Ap CT (r)=2^<A^(r), (12) 

r 

where {s°} are the singular values of B°" and where Ab° = ^2 t V t a r Ab1 and $£( r ) = 
Yliia^ia,r ( t )< i{ r ) ( t ) a{ v ) are ^ ne expansion coefficients of the change in the potential and the 
occupied-virtual orbital products transformed with the left and right singular vectors, (V° r ) 
and (U[ ar ), respectively. Here, the transformed expansion coefficients b r refer to the trans- 
formed potential basis functions g"{j) = Et^rftW- Thus, we notice that if one of these 
transformed potential basis functions g° (r) corresponds to a singular value s r that is zero 
or very small, the corresponding expansion coefficient b a r (r) can be changed (almost) freely 
without affecting the density. Therefore, an additional criterion is necessary for singling 
out the optimal optimized potential among those yielding the same density within the finite 
basis set. 



9 



2.3 Choosing the Optimal Optimized Potential 



One possibility for unambiguously singling out an optimized potential was suggested in 
Ref. [46J. This scheme starts from the requirement that the optimized potential obtained 
with a finite orbital basis set should be as close as possible to the one obtained in the basis 
set limit. Specifically, the density calculated from the optimal optimized potential should 
still agree with the target density Pq (r) or p^ (r). Thus, we introduce a complete set of virtual 
orbitals (see Ref. jl6] for details), 

OCC CT 

C(r') = 5(r-r')-^0J(r')0J(r), (13) 

i 

with the Dirac delta function S(r — r') and where the second term ensures the orthogonality 
of 0r ( r ') an d the occupied orbitals <^f (r). With this complete representation of the virtual 
orbital space, the change in the electron density due to a variation in the potential [Eq. (fTTl) ] 
can be approximated as 

OCC CT 

Ap°(r) «X>f(r)(£|f + <|#>. ( 14 ) 

i 

For singling out the optimal potential, we search for the potential for which the electron 
density does not change considerably when the orbital basis is enlarged, and minimize 

* / \ 9 p i r occ<T i 2 

d 3 r . f 1 



£>f(r)(#|T + <|#> 



d r — > min . (15) 



p(r) ' J p a (r) 

Here, the inverse density has been introduced as a weighting function (i.e., the relative change 
in the density is minimized) to obtain a uniformly accurate potential. As is discussed in 
Ref. [16], this choice can also be justified using theoretical arguments. The minimization 
then leads to the linear systems of equations A cr Ab cr = —z a with 



and 
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where hZ 



A/2 + Vg{r). This problem can be solved directly, without explicitly solving 



the KS equations using an extended orbital basis set. Results obtained using this scheme 
will be referred to as 11 optimal (full)" in the following. 

For comparison, we will also employ two additional schemes for singling out one optimized 
potential. The first one, called "balanced" in the following, is based on the idea that a unique 
potential is also obtained if the potential basis set is chosen such that it is balanced with 
respect to the orbital basis set [22] • This can be achieved by only retaining those transformed 
potential basis functions g°(r) corresponding to singular values s r that are not too small, 
i.e., above a chosen threshold s t hr- This is closely related to the OEP scheme of Kollmar 
and Filatov [23]. Note that in the spin- unrestricted case considered here, such a scheme 
effectively employs different potential basis sets for the a- and /3-spin potentials. 

In addition, we also use a criterion for singling out the optimized potentials that are as 
smooth as possible (labelled " smooth" in the following). To this end, we minimize the norm 
of the gradient of the potential, 



under the constraint that the change in the density [Eq. (|12|) ] is below a chosen threshold 
e t h r . This results in a quadratic programming problem that can be solved using standard 
approaches [46J. This criterion is in close analogy to the method of Yang and coworkers 
[41,45,62], who introduced a similar constraint by employing a penalty function during the 
direct optimization. Note that a common feature of all three approaches presented here is 
that they are applied a posteriori, and hence a direct optimization of the potentials must 
be performed first. This first step then provides a non-unique potential and corresponding 
orbitals, which are required for the following second step. 



/ 



V<(r)| 2 d 3 r 




(18) 
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3 Computational Methodology 



All finite basis set calculations were performed with a local version of the Amsterdam Density 
Functional (ADF) program package [53] together with the PyAdf scripting framework [64J. 
To allow for the treatment of spin-unrestricted target densities, we extended our recent 
implementation [30] of the Wu-Yang direct optimization algorithm and of the subsequent 
step for singling out an unambiguous optimized potential [16]. 

The TZ2P and QZ4P Slater-type orbital (STO) basis sets of ADF were used as orbital basis. 
The potential was expanded in a finite basis set [Eq. (jSJ], using ADF's density fitting basis 
sets corresponding to the TZ2P or QZ4P orbital basis sets. In all calculations, these basis sets 
were augmented with additional Is functions in an even-tempered fashion (see Supporting 
Information for details). As initial guess for the potential, we used a scaled version of the 
Fermi-Amaldi potential JHJ1IS3] of the fixed target electron density Po( r ) = Po ( r ) + Po( r )) 
namely 



where £ represents the most diffuse exponent in the STO orbital basis set. This scaled 
Fermi-Amaldi potential ensures that the optimized local potentials have the correct long- 
range behavior. In the case of target densities obtained from wave-function based ab initio 
calculations, which employed Gaussian-type orbital (GTO) basis sets, the Coulomb potential 
^Coui[po] and Fermi-Amaldi potential t>FA[po] in Eq. ([8]) are evaluated for an approximate 
reference density po(r) obtained from a DFT calculation in ADF in which the orbitals are 
expanded in STO basis functions. 

If large basis sets are employed for the potential, the Hessian matrix of Eq. ([TO]) contains many 
small eigenvalues which decay gradually to zero. This causes convergence problems during 
the Newton-Raphson optimization, which we previously addressed by ignoring eigenvalues 
below a certain threshold. However, for ab initio target densities expanded in GTOs, this 
scheme still caused poor convergence behavior. Therefore, we followed the work of Wu and 




(19) 
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Yang [05] and performed a singular value decomposition (SVD) of the Hessian H. Then the 
inverse Hessian can be expressed as 

H 1 = U diag(l/ oy) V T , (20) 

where the columns of U and V are the left and right singular vectors, respectively, for the 
corresponding singular values s r . To this inverse Hessian, a Tikhonov regularization [5^157] 
is applied by replacing it by 

H- 1 = Udiag(/ r /(7 r )V T , (21) 

where f r is a filter factor, which is chosen as 

We found that appropriate values for the parameter A turn out to be 10 -4 < A < 0.01. If 
<y T ^> A, the filter factor f r is approximately one, while in the case of cry <C A, f r approaches 
zero. Thus, instead of abruptly discarding small singular values, the Tikhonov regularization 
cuts them off gradually. 

As convergence criterion for the Wu-Yang direct optimization, we used the absolute error 
in the a- and /3-electron densities A abs compared to the target a- and /3-electron densities, 
defined as 

A: bs = / |p CT (r)-/^(r)|d 3 r. (23) 

The minimal absolute error that can be achieved depends on the considered system, as will 
be discussed below. 



When singling out the optimal potential according to the scheme of Ref. [46J, regions in 
which the electron density is very small turn out to be problematic. This is because for 
densities expanded in finite GTO or STO basis sets, the exact potential in these regions 
shows artifacts caused by unphysical nodes in the density [68,69] . To avoid these artifacts, 
grid points at which the reconstructed a- or /3-electron densities are smaller than a threshold 
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are ignored when constructing the right-hand side z a according to Eq. (|17|) . This corresponds 
to assuming that the optimized potentials are already well approximated by the initial guess 
if the difference between the target and reconstructed density is smaller than the threshold. 
For the lithium atom, this threshold was chosen as 10~ 4 , whereas for the dioxygen molecule 
it was set to 1CT 8 . 

For target densities expanded in STOs, the exchange-correlation potential is obtained by 
adding the Fermi-Amaldi potential to the part of the potential expanded in basis functions 
[cf. Eq. OH])]. In the case of target densities obtained in GTOs, the final exchange-correlation 
potentials for a- and /3-electrons are obtained as 

< c (r) = (^CouiNO) - vcoui[pw Y ](r)) + v FA [p Q ](r) + £ b a t g t {v). (24) 

t 

Here, the first term accounts for the difference between the Coulomb potential used as initial 
guess (evaluated for the density p ) and the Coulomb potential corresponding to the target 
density p . The latter is approximated by the density pwy obtained in a STO expansion 
from the Wu-Yang optimization before singling out an unambiguous potential, as this density 
represents to the best available STO representation of the GTO target density. 

To obtain numerical reference potential for atoms, we employed a modified van Leeuwen- 
Baerends algorithm [16] in combination with a numerical solution of the KS equations on 
a logarithmic radial grid [70l[7T], as described in Ref. [16]. Here, we used the same initial 
guess for the potential, and updated the potential iteratively until the absolute error A^ bs 
compared to the target a- and /3-electron densities was each below 10~ 4 ebohr -3 . 

All CASSCF calculations for obtaining accurate ab initio target densities were performed 
with the Molpro program package [72] using Dunning's cc-pVTZ basis set for all atoms 
[?3J[7I] . For the lithium atom, all electrons are correlated in all orbitals [corresponding to a 
full configuration interaction (Full-CI) treatment], while for the oxygen molecule the electron 
(spin) density from a CAS(12,12)SCF calculation was employed. Here, we verified that the 
resulting densities are converged with respect to the dimension of the active space. 
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4 Optimized Potentials from Spin Densities 



4.1 The Lithium Atom 

BP86 target (spin) density expanded in STOs 

As a simple test case, we consider the lithium atom. In its doublet ground-state, there are 
two a-electrons and one /3-electron, i.e., one unpaired electron. First, we use the a- and j3- 
electron densities from a unrestricted KS-DFT calculation employing the QZ4P orbital basis 
set and the BP86 exchange-correlation functional as target. Here, it should be possible to 
reconstruct the target (spin) density accurately in a potential reconstruction using the same 
orbital basis set. The target total and spin densities are shown in Fig. [TJ For the lithium 
atom, there is only minimal spin polarization and the spin density is determined by the 
unpaired electron in the 2s orbital, whereas almost identical Is orbitals are obtained for a- 
and /3-electrons (see also Fig. 1 in the Supporting Information). 

[Figure 1 about here.] 

For assessing the quality of the optimized potentials obtained with finite orbital basis sets, 
we determined the exchange-correlation potentials for a- and /3-electrons numerically as 
described above. These are presented in Fig. [2] as reference (black line) alongside the poten- 
tials obtained from potential reconstruction calculations employing the finite QZ4P orbital 
basis set and the corresponding QZ4P density fitting basis set augmented with additional 
tight Is functions for expanding the potential. In addition, we included the BP86 exchange- 
correlation potential evaluated for the target density (blue dashed line), i.e., the potential 
that was used in the finite-basis set KS-DFT calculation for determining the target a- and 
/3-electron densities. We note that, even though it is close to it, this BP86 potential does not 
agree with the numerical reference potential. As was pointed out before, these two potentials 
should only be equal in the basis set limit [46|l69]. 
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[Figure 2 about here.] 



The reconstructed a- and /3-electron exchange-correlation potentials v° c {v) and v@ c (r) are 
shown in Figs. [2^ and [2b. In both cases, the potentials obtained from the Wu-Yang opti- 
mization in the first step show large oscillations and are, therefore, not shown in the figures. 
These oscillations are removed if an unambiguous potential is singled out in the second step. 
Irrespective of which of the schemes described in Section [2T31 is applied, the potentials closely 
agree with the numerical reference for r > 0.5 bohr. However, differences are found closer to 
the nucleus (see also the insets in Fig. [2]). When using an implicitly balanced potential basis 
set (dashed red line), the potential is too large close to the nucleus. The a-electron potential 
has a spurious minimum in this case and show slight oscillations around the numerical refer- 
ence potential. Singling out the potential that is as smooth as possible (solid magenta line) 
does not introduce oscillations, but also results in a potential with the wrong behavior for 
small r. On the other hand, the optimal potential determined using the criterion of Eq. f fT5|) 
closely matches the numerical reference potential also close to the nucleus. 



[Table 1 about here.] 



For a more quantitative comparison of the different approaches, we present the absolute errors 
A^" um [cf. Eq. (l2U|) ] in the a- and /3-electron densities obtained from the different potentials in 
a numerical solution of the KS equations compared to the target a- and /3-electron densities 
in Table HI Naturally, this absolute error is the smallest for the numerical reference potential. 
For the a-electron potentials, the largest error is obtained for the potential obtained from the 
Wu-Yang optimization with ca. 10~ 2 e bohr -3 , i.e., the error is two magnitudes larger than for 
the numerical reference. This error is only slightly reduced by using an implicitly balanced 
basis set or by singling out a smooth potential. In contrast, for the optimal potential, 
the absolute error in the numerical density is reduced by an order of magnitude below 
10~ 3 ebohr -3 . A similar picture is obtained for the /3-electron potentials, even though all 



16 



the errors are smaller in this case. Thus, these results confirm the previous finding that high- 
quality potentials can be obtained by applying the criterion of Eq. f JT5|) for unambiguously 
singling out the optimal potential [46J . 

In addition to the errors in the densities obtained from a numerical solution of the KS 
equations with the different potentials, Table [I] also includes the absolute errors A^ mte 
obtained with these potentials in the finite orbitals basis set. In this case, the smallest error 
is obtained for the potentials obtained directly from the Wu-Yang optimization, and these 
absolute errors correspond to the convergence criterion used in this step. After singling 
out one optimized potential in the second step the error increases, but the smallest one 
is obtained for the optimal potential. Note that for both the optimal and the numerical 
reference potential the absolute density errors are larger in the finite orbital basis set than for 
the numerical solution of the KS equations. This discrepancy was discussed previously [IB] 
and arises because it is in general not possible to reproduce the target density both in a 
given finite basis set and in a fully numerical calculation at the same time. However, as was 
noted before, these differences decrease when increasing the size of the orbital basis set. 

After assessing the quality of the reconstructed a- and /3-electron exchange-correlation po- 
tentials, we turn to Figs. |2b and Figs. EH, where the same results are shown as total exchange- 
correlation potentials v^(r) = |(t>x C ( r ) + v L( r )) an d spin exchange-correlation potentials 
^xc 111 ! 1 ") = \( v xc( r ) ~~ v xc( r ))i respectively. While for the total potential, the overall results 
are similar to those discussed for the individual spin components, it is apparent that re- 
constructing the spin potential accurately is significantly more difficult. To calculate the 
difference between the a-electron potential and the /3-electron potential reliably, it is nec- 
essary to determine each of these individual a- and /3-electron potentials with comparable 
accuracy. 

For the lithium atom considered here, the spin polarization is very small. Therefore, the spin 
part of the exchange-correlation potential is very small as well. In particular, it is almost 
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constant in the range probed by the Is orbital and, therefore, does not introduce a significant 
spin polarization for the Is orbital. Thus, the spin potential is mostly due to the different 
asymptotic decay of the a- and /3-electron densities. Since the BP86 exchange-correlation 
potentials evaluated from the target a- and /3-electron densities do not show the correct 
asymptotic decay, the corresponding spin potential almost vanishes. Even though the spin 
potential is evaluated as the difference of the much larger a- and /3-electron potentials, the 
optimal potential reproduces the numerical reference almost perfectly. On the other hand, 
the smooth and the balanced potential deviate from the numerical reference not only close to 
the nuclei (where these differences were also recognizable for the a- and /3-spin potentials), 
but also further away from the nuclei. Thus, in order to reconstruct the spin potential 
t>^P m (r) reliably, it is essential to single out the optimal potential according to the criterion 
of Eq. (1T5|). 

Full-CI target (spin) density expanded in GTOs 

Next, we use the accurate a- and /3-electron densities from a Full-CI calculation for the 
lithium atom as our target. In this case, the target a- and /3-electron densities have to 
be expanded in a GTO orbital basis set, which might result in additional difficulties when 
performing the potential reconstruction with our implementation using an STO orbital basis 
set. The Full-CI total and spin densities are included in Fig. [T] and are on the scale of 
the plots almost indistinguishable from the BP86 ones considered above. Nevertheless, the 
difference densities also included in Fig. [TJ show that there are slight differences, in particular 
for the Is orbital density and the asymptotic decay. 

[Figure 3 about here.] 

Fig. |3] shows the reconstructed total and spin exchange-correlation potentials, whereas the 
individual a- and /3-electron potentials are presented in Fig. 2 in the Supporting Information. 
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The potentials obtained from the Full-CI density with a numerical solution of the KS equa- 
tions are included as black dashed lines in these plots. These reference potentials features 
considerable oscillations which are most pronounced near the nucleus. Such an oscillatory 
behavior in the reconstructed potentials is commonly found when the target density is ex- 
panded in a set of Gaussian-type orbitals (GTOs) [75] and can be attributed to the deficiency 
of GTOs to represent the correct form of the electron density close to the nucleus. However, 
these oscillations can be reduced when the GTO orbital basis set used for determining the 
target density is enlarged [75J. This is shown in Fig. 3 in the Supporting Information by 
comparing numerical potentials reconstructed from the Full-CI/cc-pVTZ target density and 
from a CAS (3,30) SCF/cc-pVQZ target density. 

As a first step of the potential reconstruction, we determine a (non-unique) potential using 
the Wu-Yang direct optimization algorithm. Here, we apply two different STO orbital 
basis sets (TZ2P and QZ4P) in combination with the corresponding density fitting basis 
sets augmented with additional tight Is functions for expanding the potential. While for 
the BP86 target density expanded in STOs this optimization could be converged until the 
absolute error in the density was smaller than 8.0 ■ 10" 4 and 10~ 4 ebohr -3 for the a- and 
/3-electron density, respectively, looser convergence criteria have to be used for the GTO 
target densities. With a TZ2P orbital basis set, the absolute density error in the Wu- 
Yang optimization only reaches 5.1 ■ 10 -3 and 2.0 ■ 10 -3 ebohr -3 for the a- and /3-electron 
density, respectively (see Tab. HT|) . This is further reduced to 2.5 • 10 -3 and 4.0 • 10~ 4 e bohr -3 , 
respectively, with the larger QZ4P orbital basis set. The corresponding difference densities 
are shown in the upper insets in Fig. [TJ Thus, a sufficiently large STO orbital basis set is 
required to be able to represent the target density from an ah initio calculation using GTOs. 
Nevertheless, with the QZ4P basis set it is possible to achieve an agreement close to the 
one obtained for the target density expanded in STOs. We note that the methodological 
improvements discussed in Section [3j in particular the use of the Tikhonov regularization, 
are essential here to make this convergence possible. 
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[Table 2 about here.] 



Since the numerical potential reconstructed from the GTO densities show oscillations due to 
the insufficiencies of the GTO basis set close to the nucleus, we also performed a numerical 
potential reconstruction using the a- and /3-electron densities from the Wu-Yang optimiza- 
tion, which are expanded in an STO orbital basis, as target. These are included in Fig. [3] 
as solid black line. They do not show oscillations near the nucleus anymore, but otherwise 
closely match the numerical potentials reconstructed from the GTO densities. However, for 
the TZ2P orbital basis set, the total and spin potentials reconstructed from the STO and 
GTO densities differ significantly for r > 2 bohr. This can be traced back to the /3-electron 
potential (see Fig. 2 in the Supporting Information), and is caused by spurious nodes ap- 
pearing in the /3-electron density in the region where it is very small [68,69j. The /3-electron 
density in the lithium atom is particularly prone to such artifacts because it is due to only 
one orbital, i.e., nodes in this orbital are not offset by other (nonzero) orbitals. In the fol- 
lowing, we will consider the numerical potentials reconstructed from the Wu-Yang densities 
expanded in STOs as reference for the finite-basis set potential reconstruction. 

The optimized potentials obtained from the finite-basis set reconstruction using the different 
schemes for singling out one unambiguous potential after the Wu-Yang optimization are 
included in Fig. [3j The potentials obtained directly from the Wu-Yang optimization contain 
considerable oscillations and are hence not shown in the figure. First, we consider the total 
exchange-correlation potentials in Figs. [3^l and |3fc (for the individual a- and /3-electron 
potentials, see Fig. 2 in the Supporting Information). With the TZ2P orbital basis set, both 
the smooth and the optimal potential show a good agreement with the numerical reference 
potential for r < 2 bohr, whereas the potential obtained with an implicitly balanced potential 
basis set still features some small oscillations. For r < 2 bohr, none of the potentials 
can reproduce the bump in the reference potential caused by artifacts of the /3-electron 
density in this region. In principle, the optimal potential should account for this feature, 
since the criterion of Eq. ffT5]) minimizes the deviation from the Wu-Yang density in the 
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basis-set limit, i.e., it should converge to the numerical reference potential determined from 
that density. However, as discussed in Section |3l regions where the density is smaller than 
10 -4 ebohr -3 , which is reached for the /3-electron density at r > 2 bohr, are discarded in 
our implementation to avoid such artifacts at small densities. In fact, we verified that after 
increasing the threshold for discarding small densities, the optimal potential does account 
for the bump in the reference potential. 

When increasing the size of the orbital basis set to QZ4P, the total exchange-correlation 
potentials from all three schemes closely agree with the numerical reference. Differences 
only remain for the smooth and balanced potentials close to the nucleus, while the optimal 
potential matches the reference also in this region. For a more quantitative comparison, 
the absolute errors A^° um in the a- and /3-electron densities obtained from the different 
potentials in a numerical solution of the KS equations compared to the target a- and (3- 
electron densities are listed in Table [TTJ In all cases, the smallest error is achieved for the 
optimal potential. Moreover, the absolute errors in the numerical densities decrease for all 
three schemes when going from the TZ2P to the QZ4P orbital basis set. 

In addition, Table [TT] includes the absolute errors A^ nite in the a- and /3-electron densities 
obtained with the different potential in the finite orbital basis set. In general, these errors 
increase compared to the potentials obtained from the Wu-Yang procedure when applying 
the schemes for singling out one optimized potential. For the TZ2P orbital basis, the largest 
error in the a-electron density within the finite basis set is in fact obtained for the optimal 
potential, even though for this potential the error in the numerical density is the smallest. 
Again, this is because the error in the finite-basis set density and in the numerical density 
cannot be minimized at the same time [IB]. Note that the finite-basis set density error 
with the optimal potential is comparable to the one obtained with the numerical reference 
potential. When going to the larger QZ4P orbital basis set, the finite-basis set density 
errors decrease significantly, both for the optimal and for the numerical reference potentials. 
This can also be seen in the finite-basis set difference densities obtained with the optimal 
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potentials, which are included in the lower insets in Fig. [3J 

Finally, we turn to the reconstructed spin exchange-correlation potentials shown in Fig. [3b 
and [3H. Again, reproducing the numerical reference potential is much more difficult in this 
case because the spin potential is calculated as the difference of the a- and /3-electron poten- 
tials. Therefore, even though with the TZ2P orbital basis set, the smooth potential visually 
agrees with the reference for the total and the individual a- and /3-electron potentials at 
r < 2 bohr (see Fig. [3^ and Figs. 2a and 2b in the Supporting Information), the correspond- 
ing spin potential deviates significantly. On the other hand, the optimal potential shows a 
good agreement with the numerical reference for r < 2 bohr and only fails to reproduce the 
spurious behavior for larger r (see discussion above). With the QZ4P orbital basis set, the 
agreement with the numerical reference improves for all reconstructed potentials. However, 
for the potential obtained with an implicitly balanced basis set, there are still significant 
oscillations, which are not obvious in the plots of the individual a- and /3-electron potentials 
(see Figs. 2c and 2d in the Supporting Information). The smooth potential qualitatively 
reproduces all features of the numerical reference, but the best agreement is achieved for the 
optimal potential. 

The comparison of the TZ2P and the QZ4P orbital basis sets demonstrates that the optimal 
potentials obtained with the smaller TZ2P orbital basis already agree very well with the 
numerical reference potential determined from the Wu-Yang density, which corresponds to 
the closest available approximation of the target density in the STO orbital basis. However, 
this Wu-Yang density differs from the target density expanded in GTOs so that deviations to 
the numerical potential determined from that target density occur. Thus, the QZ4P orbital 
basis set is required in order to obtain a (spin) potential that agrees with the numerical 
potential calculated from the GTO target density not because of the basis set requirements 
of the scheme for singling out the optimal potential, but because of the need to reproduce 
the target density represented in an GTO basis set with STOs in the Wu-Yang optimization 
step. However, the use of an STO representation of the density in the first step before 
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determining the optimal potential according to the criterion of Eq. ( I15p has the advantage of 
avoiding the spurious oscillations in the reconstructed potential arising for a GTO expansion 
of the target density close to the nucleus. Of course, a sufficiently large STO basis set will 
be able to reproduce the spurious behavior of the GTO density close to the nucleus and will 
result in a potential showing the corresponding oscillatory behavior. 

Finally, we note that the optimal potentials presented here are converged with respect to 
the size of the potential basis set. Adding additional tight or diffuse functions does not 
alter the resulting optimal potentials significantly anymore. The dependence of the optimal 
potentials on the potential basis set is illustrated in Figs. 4 and 5 and Table I and II in the 
Supporting Information. 

4.2 An Open-Shell Molecule: Dioxygen 
BP86 target (spin) density expanded in STOs 

We now investigate a diatomic molecule with an open-shell ground state, namely dioxygen O2 
with an 0-0 bond distance of 1.21 A in its equilibrium structure. Here, the two antibonding 
7r*-orbitals are singly occupied. First, we consider a target density from an unrestricted KS- 
DFT calculation using the BP86 exchange-correlation functional and the Slater-type QZ4P 
orbital basis set. This target spin density is shown in Fig. and [lb along the bonding axis 
( x-axis) and perpendicular to it (y-axis), respectively. Along the y-axis, the spin density is 
determined by the singly occupied orbitals which have a cylindrical shape around the bond 
axis, as is also visible in the plot of the spin density in the xy-plane in Fig. Hfc. On the other 
hand, the singly occupied orbitals vanish on the bond axis and, therefore, the spin density 
along the 2- axis is solely due to the spin polarization of the doubly occupied orbitals (i.e., 
the differences between the occupied a- and /9-electron orbitals). In particular, there is a 
region where the spin density becomes negative (see the yellow regions in Fig. Hb). 
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[Figure 4 about here.] 



The total exchange-correlation potentials reconstructed for the dioxygen molecule from the 
BP86/QZ4P target density using the different scheme for singling out an unambiguous po- 
tential are shown in Fig. and b along the bond axis and perpendicular to it, respectively. 
For comparison, the figures include the BP86 exchange-correlation potential calculated from 
the target density, i.e., the exchange-correlation potential used for determining the target 
density. Note, however, that because the finite QZ4P basis set is employed when calculating 
the target density, this BP86 potential is not equal to the exact potential corresponding to 
the target density [33ll4"6] . 

[Figure 5 about here.] 

The different reconstructed total exchange-correlation potentials all agree rather accurately 
with the BP86 exchange-correlation potential on the scale of the plots. Recognizable differ- 
ences are only observed along the bond axis in the region between the oxygen atoms. The 
BP86 potential shows a plateau in this bond region, which is also reproduced when singling 
out the optimized potential that is as smooth as possible. On the other hand, the balanced 
as well as the optimal potential do not exhibit such a plateau, but have a maximum at 
the midbond point. Similar observations can be made for the individual a- and /3-electron 
potentials (see Fig. 7 in the Supporting Information). 

[Table 3 about here.] 

The absolute errors A^ mte in the a- and /3-electron densities obtained from the different 
reconstructed potentials within the finite orbital basis set compared to the target densities 
are listed in Table IIHI Since the QZ4P orbital basis set used for representing the target 
density is also employed in the potential reconstruction, the Wu-Yang optimization can be 
converged such that the error in the a- and /3-electron densities is below 2.0 • 10 -4 ebohr~ 3 . 
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This error increases by one order of magnitude for the potential obtained with an implicitly 
balanced potential basis set and further increases when determining the optimized potential 
that is as smooth as possible. Note that these errors depend on the choice of the threshold 
for discarding small singular values s t h r and for the change in the density e t hr, respectively. 
For the optimal potential, the finite-basis set absolute error in the a- and /3-electron densities 
increases further to 2.8 ■ 10~ 2 and 2.3 ■ 10~ 2 ebohr -3 , respectively. The differences between 
the reconstructed total and the spin densities and the respective target (spin) density is 
illustrated in Fig. 14 in the Supporting Information. This comparison shows that except for 
the region close to the nuclei, all reconstructed potentials reproduce all qualitative features 
of the target spin density. 

However, to judge the quality of the different potentials, it would be necessary to determine 
the absolute errors in the a- and /3-electron densities obtained from these potentials in a 
numerical solution of the KS equations. Unfortunately, this is not easily possible for the 
molecular system considered here. Nevertheless, for the lithium atom considered above and 
the atomic systems investigated in Ref. [46J, it was demonstrated that this error is the 
smallest for the optimal potential. Thus, we expect that also for the dioxygen molecule, the 
optimal potential should be closest to the exact potential, and that the increased absolute 
errors in the finite basis set arises because it is not possible to reproduce the target density 
both in a numerical calculation and in the finite QZ4P orbital basis set. Note that the 
criterion for singling out the optimal potential [Eq. f )15p ] minimizes the error in the numerical 
densities, whereas the other approaches minimize the finite-basis set error. 

Finally, we turn to the reconstructed spin exchange-correlation potentials, which are pre- 
sented in Figs. \5jp and d. As for the lithium atom, reconstructing this spin potential is more 
difficult than reconstructing the total potential, because it is determined as the difference 
between the a- and /3-electron potentials and its accurate determination thus requires that 
these are both obtained with uniform accuracy. In the plots along the bond axis (see Fig. [5b), 
both the smooth and the balanced spin potentials deviate significantly from the BP86 spin 
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potential. In particular the latter shows rather pronounced oscillations. Note that these 
oscillations were not visible in the plots of the individual a- and /3-electron potentials, but 
they are amplified when considering the spin potential. The optimal spin potential plotted 
along the bond axis shows a more regular form and no oscillations that would appear un- 
physical, but it also differs from the BP86 spin potential. Nevertheless, except for the spikes 
at the nuclei themselves, it has a similar overall shape near the nuclei, i.e., a symmetric well 
in which the potential is smaller at the nucleus than ca. 0.5 bohr away from it. This shape 
of the spin potential corresponds to the positive spin density close to the nuclei. However, 
in the midbond region, the shape of the optimal spin potential qualitatively differs from the 
BP86 spin potential and the optimal spin potential shows a maximum at the midbond point, 
whereas the BP86 spin potential has a minimum. Here, the maximum of the optimal spin 
potential is actually in line with the negative spin density at the midbond point. Despite the 
observed differences, it appears that of the different reconstructed spin potentials, the opti- 
mal potential is closest to the BP86 spin potential. Note again that the BP86 spin potential 
differs from the exact spin potential, which is not available to us here. Furthermore, it is 
important to recall that the singly occupied orbitals have a node at the bond axis, i.e., the 
spin potential here is only responsible for the rather small spin polarization of the doubly 
occupied orbitals (see Fig. H^) and, therefore, it should be relatively small itself. 



[Figure 6 about here.] 



When inspecting the reconstructed spin potentials perpendicular to the bond axis (see 
Fig. EH), all reconstructed spin potentials qualitatively agree with the BP86 spin poten- 
tial. Here, the spin potential is significantly larger than along the bond axis, as it now covers 
the region to which the singly occupied orbitals extend (cf. Fig. 0b). While the differences 
are small, the optimal spin potential is closest to the BP86 spin potential and only differs at 
the nucleus, whereas both the smooth and the balanced spin potentials show some oscillatory 
behavior. The overall shape of the optimal and the BP86 spin potentials in the xy-plane 
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is compared in Fig. Except for the differences along the bond axis, in particular in the 
midbond region, already discussed above, we observe a good overall agreement. With the 
results obtained for the lithium atom in mind, we attribute these differences to the fact that 
the BP86 potential does not agree with the exact spin potential, and expect that the optimal 
spin potential is actually a closer approximation to the exact spin potential. 

CASSCF (spin) density expanded in GTOs 

After considering the reconstruction of the spin exchange-correlation potential for the dioxy- 
gen molecule for a BP86/QZ4P target (spin) density, we now turn to an accurate target (spin) 
density obtained from a CAS(12,12)SCF/cc-pVTZ calculation. This ab initio spin density is 
compared to the BP86 spin density in Fig. HI Perpendicular to the bond axis (see Fig. Hb), 
where the spin density is dominated by the unpaired 7r*-electrons, the CASSCF spin density 
qualitatively agrees with the one obtained with BP86, even though the magnitude at the 
maximum and the asymptotic form slightly differ. Larger differences are found along the 
bond axis (see Fig. H^), where the spin density is solely due to the spin polarization of dou- 
bly occupied orbitals. The CASSCF spin density is significantly larger close to the nuclei, 
whereas the negative spin density at the ends of the molecule is reduced. This is also visible 
in the plot of the CASSCF spin density in the xy-plane shown in Fig. SU. 

[Figure 7 about here.] 

While for the BP86/QZ4P target density the BP86 exchange-correlation (spin) potential 
calculated from the target density could be used for comparison — even though it is not 
identical to the exact (spin) potential — we now have no reference potential available. Nev- 
ertheless, we can still use the results obtained for the BP86 target density to judge whether 
the (spin) potentials reconstructed from the CASSCF target (spin) density are physically 
reasonable. The reconstructed total and spin exchange-correlation potentials obtained us- 
ing a QZ4P orbital basis set using the different scheme for singling out an unambiguous 
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optimized potential are presented in Fig. [7J The corresponding individual a- and /3-electron 
potentials are given in Fig. 8 in the Supporting Information. 

For the reconstructed total exchange-correlation potential, the optimal potential has a sim- 
ilar shape as the BP86 potential discussed above. Along the bond axis (see Fig. Wfr), it has 
a maximum at the midbond point and exhibits some shell structure at ca. 0.5 bohr from 
the two nuclei. A similar shape as in the outer region along the bond axis is found perpen- 
dicular to the bond axis for the optimal potential (see Fig. UJ). Both the smooth and the 
balanced total potentials have a similar shape, but the shell structure is more pronounced. 
For the balanced total potential, larger oscillations appear close to the nuclei. Nevertheless, 
the different reconstructed potentials are qualitatively rather similar. A different picture is 
obtained if the potential reconstruction is performed with the smaller TZ2P orbital basis set, 
as is shown in Fig. 11 in the Supporting Information. In this case, both the balanced and the 
smooth total potentials show large oscillations, even though the same thresholds as for the 
QZ4P orbital basis set are applied. However, the optimal potential is free of such unphysical 
oscillations already with the smaller orbital basis set. Moreover, the optimal total potential 
obtained with the TZ2P orbital basis set is on the scale of the figures in good agreement 
with the one reconstructed using the larger QZ4P orbital basis set. 

The absolute errors in the a- and /3-electron densities compared to the target densities 
obtained from these different reconstructed potentials in the finite TZ2P and QZ4P orbital 
basis sets are listed in Table IHIl Because we are now trying to reproduce target a- and 
/3-electron densities expanded in a GTO basis set with STOs, these errors are larger than 
for the BP86 target densities that was represented in the same orbital basis also used in the 
potential reconstruction. With the TZ2P orbital basis set, the Wu-Yang optimization can 
be converged until the absolute error is approximately 4 ■ 10~ 2 ebohr -3 for both the a- and 
/3-electron density. With the larger QZ4P orbital basis set, a slightly smaller error can be 
achieved for the a-electron density, whereas the error in the /5-electron density almost halved 
for both electron spins. With both orbital basis sets, the errors are larger for the smooth and 
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balanced potentials, with this increase being controlled by the chosen thresholds. For the 
larger QZ4P orbital basis set, we observe that the increase is smaller than with the TZ2P 
orbital basis set. 

The largest density errors in the finite orbital basis sets are found for the optimal potentials. 
For the TZ2P orbital basis set, the errors are approximately twice as large as for the Wu- 
Yang potentials. However, a significantly smaller increase is observed with the QZ4P orbital 
basis set. Here, we stress again that for judging the quality of the different potentials, it 
would be necessary to calculate the errors in the densities obtained with a numerical solution 
of the Kohn-Sham equations on the respective potentials. For the optimal potentials, these 
numerical density errors, which differ from the finite basis set density errors, should be 
minimized. With increasing size of the orbital basis set used in the potential reconstruction, 
the difference between the finite basis set and the numerical density errors should become 
smaller. Thus, the decrease of the finite basis set density errors when going from TZ2P to 
QZ4P supports our assumption that the optimal potentials should be closest to the exact 
potentials. The total and spin densities calculated from the different reconstructed potentials 
in the finite orbital basis set are compared to the target (spin) density in Fig. 15 in the 
Supporting Information. Here, it is obvious that these differences are larger than for the 
BP86 target (spin) density. In particular, the optimal potential yields a spin density that 
misses some of the qualitative features of the target density around the nuclei. However, it 
is also apparent that these deviations are reduced when increasing the orbital basis set from 
TZ2P to QZ4P. 

In order to study the convergence of the reconstructed potentials with respect to the poten- 
tial basis set, we performed additional calculations where additional diffuse functions were 
included in the potential basis set expansion. The resulting optimal reconstructed potentials 
are shown in Figs. 16-19 in the Supporting Information, both for the TZ2P and the QZ4P 
orbital basis set. These results show that the reconstructed potentials are converged with 
respect to the size of the potential basis set. Only far away from the molecule, i.e., in regions 
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where the density is very small, some differences between potentials expanded in different 
basis sets occur. In addition, we observe that the requirements on the potential basis set 
are most severe for the spin potential, where with the TZ2P orbital basis set increasing the 
potential basis set still has some effect. 

Finally, we now consider the spin exchange-correlation potential, which is shown in Figs. [Tfc 
and d along the 0-0 bond axis and perpendicular to it, respectively. Perpendicular to 
the bond axis, the optimal spin potential has a similar shape as the one obtained for the 
BP86/QZ4P target density with a minimum at ca. 0.5 bohr from the oxygen nuclei, i.e., where 
the spin density is the largest. In contrast, rather pronounced oscillations are observed for 
the smooth and balanced spin potentials. For the smooth and the balanced spin potentials, 
such large oscillations are also present in the plot along the bond axis, whereas the optimal 
potential is mostly well behaved. Qualitatively, it also resembles the optimal spin potential 
reconstructed from the BP86/QZ4P target density. Around the nuclei, it has negative wells, 
which are deeper for the CASSCF than for the BP86/QZ4P target density. This is in line 
with the larger spin density in this region obtained in the CASSCF calculations. In the 
bond region between the atoms, the optimal spin potential reconstructed from the CASSCF 
spin density is flatter than in the case of the BP86/QZ4P target density, which agrees with 
the smaller spin polarization in this region. Very close to the nuclei the reconstructed spin 
potential has some large oscillations. Most likely, these can be attributed to the deficiencies 
of the GTO target density, which can cause such oscillations in the exact potentials [75] . 

[Figure 8 about here.] 

To compare the overall performance of the different potential reconstruction schemes, Fig. [8] 
compares the optimal and the smooth spin potentials in the xy-plane. A similar comparison 
for the balanced spin potential is shown in Fig. 9 in the Supporting Information. While 
pronounced wiggling features around the position of the nuclei are found in the smooth 
spin potential, it is obvious that these can be eliminated in the whole xy-plane for the 
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optimal spin potential. Thus, even though there is no exact reference spin potential available 
for comparison, we find that the optimal spin potential is the only one that seems to be 
free of artifacts of the potential reconstruction. Moreover, the comparison with the spin 
potential reconstructed from the BP86/QZ4P target density shows that the optimal potential 
is physically reasonable. Therefore, we are confident that it is the reconstructed potential 
that is closest to the exact spin potential. 

Further support for this conclusion can be drawn from a comparison with the reconstructed 
potentials obtained for the same CASSCF target density with the smaller TZ2P orbital basis 
set, which are shown in Fig. 10 and Fig. 12 and 13 in the Supporting Information. Here, 
even larger oscillations are found for the smooth and balanced spin potentials, even though 
the same thresholds are applied. On the other hand, the optimal potentials are qualitatively 
similar to the QZ4P results. This is particularly obvious for the spin potential perpendicular 
to the bond axis, which is in good agreement with the one reconstructed with the QZ4P 
orbital basis set. Also along the bond axis the shapes of the optimal spin potentials are 
similar for the two orbital basis sets, even though the negative well near the oxygen nuclei 
is less pronounced. Thus, only for the optimal potential a consistent convergence with 
increasing size of the orbital basis set is observed. 

5 Conclusions 

In this work, we have extended the unambiguous reconstruction of the local potential yield- 
ing a given target density [46J to open-shell systems treated with an unrestricted KS-DFT 
formalism. Moreover, we have combined this reconstruction with the use of accurate target 
(spin) densities obtained from accurate wave-function based ab initio calculations, i.e., from 
Full-CI or CASSCF wave functions. This provides a route to accurate reference data for 
the spin exchange-correlation potential ujj£ m , which determines the spin density distribution 
Q(r) in unrestricted KS-DFT. 
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Reconstructing this spin exchange-correlation potential is a particularly challenging task, 
because it is given by the difference between the reconstructed a- and /3-electron potentials. 
Thus, one has to overcome the numerical inaccuracies caused by the ill-posed nature of 
the potential reconstruction problem in finite orbital basis sets, as these could otherwise 
be amplified when calculating the spin potential. As test cases, we chose the lithium atom 
and the oxygen molecule in its triplet state. For both systems, we considered target (spin) 
densities from unrestricted KS-DFT calculations as well as from Full-CI (for the lithium 
atom) and CASSCF (for the dioxygen molecule) calculations. These test cases made it 
possible to systematically assess the quality of the reconstructed spin exchange-correlation 
potentials. 

For the lithium atom, it is possible to compare the reconstructed potentials to the exact 
ones, which can be obtained from a fully numerical potential reconstruction. The compar- 
ison shows that the optimal spin potentials, determined using the scheme of Ref. [46J, can 
reproduce the fully numerical spin potential, while the spin potentials obtained by singling 
out the a- and /3-electron potentials that are as smooth as possible or with an implicitly 
balanced potential basis set show significantly larger deviations from the reference spin po- 
tential. In general, the quality of these smooth and balanced spin potentials strongly depends 
on the choice of the corresponding threshold values and can result in highly oscillating poten- 
tials if these are chosen too small. On the other hand, if these thresholds are too large, the 
resulting potentials lack all features present in the exact one. Moreover, different thresholds 
might be required for reconstructing the a- and /3-electron potentials in order to obtain these 
with similar quality, as is required for the reconstruction of the spin potential. 

For target densities obtained with GTO orbital basis sets, which are commonly used in wave- 
function based ab initio calculations, we reconstruct the potential using an STO orbital basis 
set. This has the advantage that the density obtained from the Wu-Yang optimization, which 
is the first step in all potential reconstruction schemes tested here, has the correct shape 
close to the nuclei and in the asymptotic region. Consequently, oscillations close to the 



32 



nucleus found in the fully numerical reference potential because of the wrong form the the 
GTO target density [75] are largely suppressed in the optimal reconstructed potentials. On 
the other hand, it requires the use of sufficiently large STO orbital basis sets in order to be 
able to reproduce the ab initio target (spin) density. For the lithium atom, we found that 
the a QZ4P orbital basis set is sufficient to obtain a good agreement with the target density. 
Moreover, we note that technical improvements to our implementation, in particular the use 
of a Tikhonov regularization in the Wu-Yang optimization [65] and the use of a cut-off value 
for discarding small density regions in the criterion of Eq. f JT5|) . were necessary to treat GTO 
target densities. 

For the dioxygen molecule, a direct comparison to a fully numerical reference spin potential 
is not possible. Nevertheless, for the target density obtained from an unrestricted KS-DFT 
calculation the exchange-correlation potential used for determining the target potential can 
provide some guidance, even though it differs from the exact potential because of the use of a 
finite orbital basis set. This comparison shows that the optimal reconstructed spin potential 
(i.e., the one obtained using the scheme of Ref. [16]) shows the best overall agreement. 
Also for the CASSCF target density, this scheme is the only one that provides a physically 
reasonable spin exchange-correlation potential, while the smooth and the balanced potentials 
are plagued by unphysical oscillations. The optimal potential shows such oscillations only 
very close to the nuclei, where they are probably due to deficiencies of the GTO expansion 
used for the target density, and to a smaller extent far away from the molecule. 

In summary, we believe that the potential reconstruction scheme proposed in Ref. [46J and 
extended here to open-shell systems provides the first reliable approach for reconstructing 
the spin exchange-correlation potential from accurate ab initio (spin) densities. The avail- 
ability of such accurate reference spin potentials can facilitate the development of improved 
spin-dependent exchange-correlation density functionals which apart from yielding accurate 
total electron densities also provide reliable spin densities. Thus, this work represents a 
prerequisite for the design of approximate exchange-correlation functional with an improved 
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spin-density dependance. 
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List of Figures 



1 Target radial (a) total densities and (b) spin densities for the lithium atom 
obtained from BP86/QZ4P and Full-CI/cc-pVTZ calculations. The difference 
between these two target densities is shown in the lower part. The insets 
present the difference between the target total and spin densities and the 
corresponding density obtained in the finite orbital basis set from the Wu- 
Yang optimized potentials (upper insets) and from the optimal optimized 
potentials (lower insets). For the Full-CI/cc-pVTZ results obtained with both 
the TZ2P and QZ4P orbital basis sets in the potential reconstruction are 
included. The corresponding plots of the a- and /3-electron densities are given 
in Fig. 1 in the Supporting Information 

2 Reconstructed potentials determined for the Li atom and a BP86 / QZ4P target 
(spin) density. The upper part shows the exchange-correlation potentials for 
(a) a electrons v" c and (b) (3 electrons v^ c , while the lower part shows (c) the 
total exchange-correlation potential v *°* and (d) the spin exchange-correlation 
potential v^. m . For the potential reconstruction with the finite QZ4P orbital 
basis set, the potentials obtained with the different schemes for singling out 
an unambiguous potential (see text for details) are shown. The accurate 
potentials obtained with a numerical solution of the KS equations ( "numerical 
(STO)") as well as the BP86 exchange-correlation potential calculated from 
the reference density ( "BP 86 xc potential") are shown for comparison. The 
latter is shifted such that it agrees with the numerical reference at r = 5 bohr 
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Reconstructed potentials determined for the Li atom and a Full-CI/cc-pVTZ 
target (spin) density with a TZ2P and with a QZ4P orbital basis set in 
the potential reconstruction. The upper part shows the (a) total exchange- 
correlation potentials and (b) the spin exchange-correlation potentials 
t>xc in obtained with the TZ2P orbital basis set in the potential reconstruction, 
whereas the lower part (c) and (d) shows the corresponding results for the 
QZ4P orbital basis set. The accurate potentials obtained with a numerical 
solution of the KS equations ( "numerical (GTO)") from the GTO target den- 
sity as well as from the Wu-Yang reconstructed density ( "numerical (STO)") 
are shown for comparison. These reference potentials have been shifted such 
that they agree with the optimal potential at r = 2.7 bohr. The corresponding 
plots of the a- and /3-electron potentials are given in Fig. 2 in the Supporting 
Information 

Target spin densities for the dioxygen molecule obtained from BP86/QZ4P 
and CAS(12,12)SCF/cc-pVTZ calculations. Both spin densities are compared 
(a) along the bond axis (x-axis), (b) perpendicular to the bond axis (y-axis). 
Furthermore, (c) and (d) show the BP86/QZ4P and CAS(12,12)SCF/cc- 
pVTZ spin densities, respectively, in the xy-plane. The corresponding plots 
of the total and the individual a- and /3-electron densities are given in Fig. 6 
in the Supporting Information 
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Reconstructed potentials determined for the dioxygen molecule and a BP86/QZ4P 
target (spin) density. The upper part shows the total exchange-correlation po- 
tential (a) along the bond axis (x axis) and (b) perpendicular to the bond 
axis along the y axis. The lower part shows the spin exchange-correlation po- 
tential t>xc in along the (c) x axis and (d) y axis. In the potential reconstruction, 
the finite QZ4P orbital basis set was employed. For comparison, the BP86 
exchange-correlation potential calculated from the reference density ( "BP86 
xc potential") is also included. This BP86 potential is shifted such that it 
agrees with the optimal potential at x — —2 bohr or at y = +2 bohr, respec- 
tively. The corresponding plots of the individual a- and /3-electron potentials 
are shown in Fig. 7 in the Supporting Information 
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(a) Reconstructed spin exchange-correlation potential v^ m determined for the 
dioxygen molecule and a BP86/QZ4P target (spin) density in the xy-p\a,ne. 
Here, only the optimal potential reconstructed within a QZ4P orbital basis 
set is included. For comparison, (b) shows the BP86 exchange-correlation 
potential calculated from the target density 
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Reconstructed potentials determined for the dioxygen molecule and a CAS(12,12)SCF/cc- 
pVTZ target (spin) density. The upper part shows the total exchange-correlation 
potential i>*£* (a) along the bond axis (x axis) and (b) perpendicular to 
the bond axis along the y axis. The lower part shows the spin exchange- 
correlation potential v^ m along the (c) and (d) y axis. In the potential 
reconstruction, the finite QZ4P orbital basis set was employed. The corre- 
sponding plots of the individual a- and /3-electron potentials as well as the 
results obtained with the TZ2P orbital basis set are shown in Fig. 8 and 
Figs. 10 and 11 in the Supporting Information, respectively. 



Reconstructed spin exchange-correlation potential potentials determined for 
the dioxygen molecule and a CAS(12,12)SCF/cc-pVTZ target (spin) density. 
For the potential reconstruction with the finite QZ4P orbital basis set (a) the 
optimal spin potential and (b) the spin potential determined by requiring that 
the optimized potential is smooth are shown. The spin potential obtained by 
implicitly balancing the orbital and potential basis sets as well as the results 
obtained with the TZ2P orbital basis set are shown in Fig. 9 and Figs 12 
and 13 in the Supporting Information, respectively 



48 



49 



41 



Figure 1: Target radial (a) total densities and (b) spin densities for the lithium atom ob- 
tained from BP86/QZ4P and Full-CI/cc-pVTZ calculations. The difference between these 
two target densities is shown in the lower part. The insets present the difference between 
the target total and spin densities and the corresponding density obtained in the finite or- 
bital basis set from the Wu-Yang optimized potentials (upper insets) and from the optimal 
optimized potentials (lower insets). For the Full-CI/cc-pVTZ results obtained with both 
the TZ2P and QZ4P orbital basis sets in the potential reconstruction are included. The 
corresponding plots of the a- and /3-electron densities are given in Fig. 1 in the Supporting 
Information. 
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Figure 2: Reconstructed potentials determined for the Li atom and a BP86/QZ4P target 
(spin) density. The upper part shows the exchange-correlation potentials for (a) a electrons 
and (b) (3 electrons v@ c , while the lower part shows (c) the total exchange-correlation 
potential and (d) the spin exchange-correlation potential v^ m . For the potential recon- 
struction with the finite QZ4P orbital basis set, the potentials obtained with the different 
schemes for singling out an unambiguous potential (see text for details) are shown. The 
accurate potentials obtained with a numerical solution of the KS equations ( "numerical 
(STO)") as well as the BP86 exchange-correlation potential calculated from the reference 
density ( "BP 86 xc potential") are shown for comparison. The latter is shifted such that it 
agrees with the numerical reference at r = 5 bohr. 
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Figure 3: Reconstructed potentials determined for the Li atom and a Full-CI/cc-pVTZ target 
(spin) density with a TZ2P and with a QZ4P orbital basis set in the potential reconstruc- 
tion. The upper part shows the (a) total exchange-correlation potentials and (b) the 
spin exchange-correlation potentials v^ in obtained with the TZ2P orbital basis set in the 
potential reconstruction, whereas the lower part (c) and (d) shows the corresponding results 
for the QZ4P orbital basis set. The accurate potentials obtained with a numerical solution 
of the KS equations ("numerical (GTO)") from the GTO target density as well as from 
the Wu-Yang reconstructed density ( "numerical (STO)") are shown for comparison. These 
reference potentials have been shifted such that they agree with the optimal potential at 
r = 2.7 bohr. The corresponding plots of the a- and /9-electron potentials are given in Fig. 2 
in the Supporting Information. 
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Figure 4: Target spin densities for the dioxygen molecule obtained from BP86/QZ4P and 
CAS(12,12)SCF/cc-pVTZ calculations. Both spin densities are compared (a) along the bond 
axis ( x-axis), (b) perpendicular to the bond axis (y-axis). Furthermore, (c) and (d) show 
the BP86/QZ4P and CAS(12,12)SCF/cc-pVTZ spin densities, respectively, in the xy-p\a.ne. 
The corresponding plots of the total and the individual a- and /3-electron densities are given 
in Fig. 6 in the Supporting Information. 
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Figure 5: Reconstructed potentials determined for the dioxygen molecule and a BP86/QZ4P 
target (spin) density. The upper part shows the total exchange-correlation potential v (a) 
along the bond axis (x axis) and (b) perpendicular to the bond axis along the y axis. The 
lower part shows the spin exchange-correlation potential v^. m along the (c) x axis and (d) 
y axis. In the potential reconstruction, the finite QZ4P orbital basis set was employed. For 
comparison, the BP86 exchange-correlation potential calculated from the reference density 
( "BP86 xc potential") is also included. This BP86 potential is shifted such that it agrees with 
the optimal potential at x — —2 bohr or at y — +2 bohr, respectively. The corresponding 
plots of the individual a- and /3-electron potentials are shown in Fig. 7 in the Supporting 
Information. 
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Figure 6: (a) Reconstructed spin exchange-correlation potential determined for the 
dioxygen molecule and a BP86/QZ4P target (spin) density in the xy-plane. Here, only the 
optimal potential reconstructed within a QZ4P orbital basis set is included. For comparison, 
(b) shows the BP86 exchange-correlation potential calculated from the target density. 
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Figure 7: Reconstructed potentials determined for the dioxygen molecule and a 
CAS(12,12)SCF/cc-pVTZ target (spin) density. The upper part shows the total exchange- 
correlation potential v *°* (a) along the bond axis (x axis) and (b) perpendicular to the bond 
axis along the y axis. The lower part shows the spin exchange-correlation potential v^ m 
along the (c) x axis and (d) y axis. In the potential reconstruction, the finite QZ4P orbital 
basis set was employed. The corresponding plots of the individual a- and /3-electron poten- 
tials as well as the results obtained with the TZ2P orbital basis set are shown in Fig. 8 and 
Figs. 10 and 11 in the Supporting Information, respectively. 
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Figure 8: Reconstructed spin exchange-correlation potential potentials determined for the 
dioxygen molecule and a CAS(12,12)SCF/cc-pVTZ target (spin) density. For the potential 
reconstruction with the finite QZ4P orbital basis set (a) the optimal spin potential and (b) 
the spin potential determined by requiring that the optimized potential is smooth are shown. 
The spin potential obtained by implicitly balancing the orbital and potential basis sets as 
well as the results obtained with the TZ2P orbital basis set are shown in Fig. 9 and Figs 12 
and 13 in the Supporting Information, respectively. 
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Table I: Absolute errors A^ bs in the a- and /9-electron densities with respect to the target 
a- and /3-electron densities (in ebohr~ 3 ) obtained with different reconstructed potentials for 
the Li atom with BP86/QZ4P target densities. The QZ4P orbital basis set is used in the 
finite-basis set potential reconstruction. A^ mtc refers to the error in the density obtained 
from the respective potentials in the finite orbital basis set, whereas A^ um is the error for 
the density obtained from a numerical solution of the KS equations. 





a o,num 
abs 


a a, finite 
abs 


a /3,num 
abs 


a /?, finite 
abs 


Numerical (STO) 


< 10- 4 


0.0007 


< to- 4 


0.0001 


Wu-Yang 


0.0098 


0.0008 


0.0036 


< io- 4 


Balanced (sthr = 10 -2 ) 


0.0090 


0.0024 


0.0016 


0.0006 


Smooth (ethr = 10~ 2 ) 


0.0087 


0.0086 


0.0082 


0.0082 


Optimal (full) 


0.0009 


0.0014 


< io- 4 


< io- 4 
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Table II: Absolute errors A^ bs in the a- and /3-electron densities with respect to the target 
a- and /3-electron densities (in ebohr~ 3 ) obtained with different reconstructed potentials for 
the Li atom and Full-CI/cc-pVTZ target densities. Results obtained both with the TZ2P 
and with the QZ4P orbital basis set in the potential reconstruction are shown. A^ mtc refers 
to the error in the density obtained from the respective potentials in the finite orbital basis 
set, whereas A^ um is the error for the density obtained from a numerical solution of the KS 
equations. 



TZ2P QZ4P 



Full-CI/cc-pVTZ 


a a,num 
abs 


a a, finite 
abs 


a /3,num 
abs 


a /?, finite 
abs 


a a,num 
abs 


a a, finite 
abs 


a /3,num 
abs 


a /?, finite 
abs 


Numerical (GTO) 


< 1(T 4 


0.0258 


< io- 4 


0.0086 


< 10~ 4 


0.0065 


< 10" 4 


0.0010 


Numerical (STO) 


0.0051 


0.0158 


0.0020 


0.0014 


0.0025 


0.0051 


0.0005 


0.0005 


Wu-Yang 


0.0782 


0.0051 


0.0301 


0.0020 


0.0165 


0.0025 


0.0021 


0.0004 


Balanced (s t hr = 10~ 2 ) 


0.0336 


0.0079 


0.0310 


0.0020 


0.0092 


0.0074 


0.0030 


0.0004 


Smooth (e t hr = 10~ 2 ) 


0.0416 


0.0105 


0.0039 


0.0041 


0.0120 


0.0071 


0.0062 


0.0061 


Optimal (full) 


0.0220 


0.0284 


0.0021 


0.0016 


0.0061 


0.0067 


0.0005 


0.0006 
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Table III: Absolute errors A^ s m e in the a- and /3-electron densities in the finite orbital 
basis set with respect to the target a- and /3-electron densities (in ebohr~ 3 ) obtained with 
different reconstructed potentials for the dioxygen molecule. Results are shown for both the 
target densities from BP86 and from a CAS(12,12)SCF calculation and in the latter case 
using both the TZ2P and with the QZ4P orbital basis set in the potential reconstruction. 

BP86/QZ4P CASSCF/TZ2P CASSCF/QZ4P 

a a, finite a 0, finite a a, finite a 0, finite a a, finite a 0, finite 
abs abs abs abs abs abs 

0.0002 0.0002 0.0419 0.0439 0.0412 0.0268 

0.0019 0.0010 0.0521 0.0537 0.0414 0.0278 

0.0077 0.0078 0.0426 0.0462 0.0416 0.0277 

0.0284 0.0230 0.0976 0.0830 0.0640 0.0509 



Wu-Yang 

Balanced (sthr = 10 -2 ) 
Smooth (ethr = 10~ 2 ) 
Optimal (full) 
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